In [ ]:
import numpy as np
import pandas as pd
from tqdm import tqdm
import plotly.graph_objects as go
import plotly.express as px
import plotly.io as pio
pio.renderers.default = "notebook+pdf"

System Variables

In [ ]:
# Event codes
ARRIVAL=0
DEPARTURE=1
event_codes = {ARRIVAL: "ARRIVAL", DEPARTURE: "DEPARTURE"}

# Patient status codes
PENDING=0
ADMITTED=1
RELOCATED=2
REJECTED=3
DISCHARGED=4
status_codes = {PENDING: "PENDING", ADMITTED: "ADMITTED", RELOCATED: "RELOCATED", REJECTED: "REJECTED"}

# Simulation parameters
SIM_TIME = 31   # Simulation time in days

# Ward parameters
WARDS = dict(
    names = ['A', 'B', 'C', 'D', 'E', 'F'],
    lams = [14.5, 11.0, 8.0, 6.5, 5.0, 13.0],
    mu_invs = [2.9, 4.0, 4.5, 1.4, 3.9, 2.2],
    urgency_points = [7, 5, 2, 10, 5, 0],
)

NUM_WARDS = len(WARDS['names'])

RELOCATION_PROBABILITIES = np.array([
    [0.0, 0.05, 0.10, 0.05, 0.80, 0.00],
    [0.2, 0, 0.50, 0.15, 0.15, 0.00],
    [0.30, 0.20, 0, 0.20, 0.30, 0.00],
    [0.35, 0.30, 0.05, 0, 0.3, 0.00],
    [0.20, 0.10, 0.60 ,0.10, 0, 0.00],
    [0.20, 0.20, 0.20, 0.20, 0.20 ,0]
])

Penalty Computation

In [ ]:
def penalty(urgency_points, counts):
    return np.dot((counts[RELOCATED] + counts[REJECTED]), urgency_points)

Simulation Setup

In [ ]:
def simulation_setup(sim_time=SIM_TIME):

    # Patients
    patients = dict( ID=[], type=[], ward=[], status=[], U1=[], U2=[], burn_in=[] )

    # Precompute arrival and departure times for each ward
    events = dict( event_type=[], patient_ID=[], time=[] )
    patient_id = 0
    for i, (lam, mu_inv) in enumerate(zip(WARDS['lams'], WARDS['mu_invs'])):
        
        # Sample patients for ward
        clock = 0
        while clock < sim_time:

            # Sample arrival time
            U1, U2 = np.random.uniform(size=2)
            clock += -np.log(U1)/lam

            # Add patient to event list if arrival time is before simulation end
            if clock <= sim_time:
                # Event
                events['time'] += [clock, clock-np.log(U2)*mu_inv]
                events['patient_ID'] += [patient_id]*2
                events['event_type'] += [ARRIVAL, DEPARTURE]

                # Patient
                patients['ID'].append(patient_id)
                patients['type'].append(i)
                patients['ward'].append(i)
                patients['status'].append(PENDING)
                patients['U1'].append(U1)
                patients['U2'].append(U2)
                patients['burn_in'].append(0)

                # Update patient ID
                patient_id += 1


    # Sort events by time
    idx = np.argsort(events['time'])
    for key in events.keys():
        events[key] = [events[key][i] for i in idx]

    return patients, events

Inspect the Event List

In [ ]:
# Setup simulation
patients, events = simulation_setup()

# Create dataframes and map codes to names
patient_df = pd.DataFrame(patients); patient_df['status'] = patient_df['status'].map(status_codes)
patient_df['type'] = patient_df['type'].map({i: WARDS['names'][i] for i in range(NUM_WARDS)})
event_df = pd.DataFrame(events); event_df['event_type'] = event_df['event_type'].map(event_codes)

# Print dataframes
print("Patients:")
print(patient_df.head().to_string(index=False), "\n")
print("Events:")
print(event_df.head().to_string(), "\n")
Patients:
 ID type  ward  status       U1       U2  burn_in
  0    A     0 PENDING 0.675212 0.036522        0
  1    A     0 PENDING 0.914004 0.612931        0
  2    A     0 PENDING 0.268128 0.037591        0
  3    A     0 PENDING 0.582789 0.018506        0
  4    A     0 PENDING 0.414830 0.484123        0 

Events:
  event_type  patient_ID      time
0    ARRIVAL        1011  0.017203
1    ARRIVAL           0  0.027085
2    ARRIVAL         451  0.029163
3    ARRIVAL           1  0.033286
4    ARRIVAL         765  0.038836 

Simulation Loop

In [ ]:
def simulate(patients, events, capacities, burn_in_period=15):

    # Clear patient status
    N = len(patients['ID'])
    patients['ward'] = [None]*N
    patients['status'] = [PENDING]*N

    # List for saving ward states
    states = []

    # Dict for saving counts
    counts = {
        ADMITTED: np.zeros(NUM_WARDS, dtype=int),
        RELOCATED: np.zeros(NUM_WARDS, dtype=int),
        REJECTED: np.zeros(NUM_WARDS, dtype=int),
    }

    # Simulate loop
    occupancy = np.zeros(NUM_WARDS, dtype=int)
    for i, (event_type, patient_id) in enumerate(zip(events['event_type'], events['patient_ID'])):

        # Check if patient is in burn-in period
        if events['time'][i] < burn_in_period:
            patients['burn_in'][patient_id] = True

        # Get patient type
        patient_type = patients['type'][patient_id]
        
        if event_type == ARRIVAL:

            # Admit patient
            if occupancy[patient_type] < capacities[patient_type]:
                ward = patients['type'][patient_id]
                status = ADMITTED

            else:
                # Relocate patient (or reject if alternative ward is full)
                ward = np.random.choice(NUM_WARDS, p=RELOCATION_PROBABILITIES[patient_type])
                status = RELOCATED if occupancy[ward] < capacities[ward] else REJECTED

            # Update patient
            if status != REJECTED:
                occupancy[ward] += 1
                patients['ward'][patient_id] = ward
            patients['status'][patient_id] = status

            if events['time'][i] >= burn_in_period:
                counts[status][patient_type] += 1

        # Discharge patient
        elif patients['status'][patient_id] != REJECTED:
            occupancy[patients['ward'][patient_id]] -= 1
    
        states.append(occupancy.copy())

    events['states'] = states

    # Check if simulation ended with non-zero occupancy
    if any(occupancy != 0):
        print('Error: Simulation ended with non-zero occupancy')

    return dict(states=states, counts=counts, penalty=penalty(WARDS['urgency_points'], counts))

Run Simulation

In [ ]:
patients, events = simulation_setup()
capacities = [49, 28, 22, 16, 16, 34]
sim_out = simulate(patients, events, capacities)
states = sim_out['states']

fig = go.Figure()
for i, state in enumerate(np.array(events['states']).T):
    fig.add_trace(go.Scatter(x=events['time'], y=state, mode='lines', name=WARDS['names'][i]))
fig.update_layout(title='Ward occupancies', xaxis_title='Time', yaxis_title='Occupancy', width=800)
fig.show()

Gradient Computation

In [ ]:
def compute_gradients(capacities):
    M = len(capacities)
    grads = np.zeros(M)
    n = 10

    def simulate_penalty(capacity_adjustment):
        total_penalty = 0
        for _ in range(n):
            capacities[i] += capacity_adjustment
            patients, events = simulation_setup()
            sim_out = simulate(patients, events, capacities)
            total_penalty += sim_out['penalty']
            capacities[i] -= capacity_adjustment
        return total_penalty / n

    for i in range(M):
        p_small = simulate_penalty(-1)
        p_large = simulate_penalty(1)
        grads[i] = (p_large - p_small) / 2

    return grads


def gradient_descent(capacities, beds_to_F=34):

    penalties = np.zeros(beds_to_F+1)
    F_rel_prob = np.zeros(beds_to_F+1)

    # Use the samee simulations to estimate performance metrics
    sim_setup_samples = [simulation_setup() for _ in range(100)]

    # Realocate beds to F
    for i in tqdm(range(beds_to_F+1)):
        
        # Simulate and compute performance metrics
        sim_out = [simulate(patients, events, capacities) for patients, events in sim_setup_samples]
        penalties[i] = np.mean(np.array([out['penalty'] for out in sim_out]))
        num = [out['counts'][RELOCATED][-1] + out['counts'][REJECTED][-1] for out in sim_out]
        den = [out['counts'][ADMITTED][-1] + out['counts'][RELOCATED][-1] + out['counts'][REJECTED][-1] for out in sim_out]
        F_rel_prob[i] = np.mean([num/den for num, den in zip(num, den)])

        # Compute gradients and update capacities
        grads = compute_gradients(capacities)
        idx = np.argmax(grads[:-1])
        capacities[idx] -= 1
        capacities[-1] += 1

    # Correct for last iteration
    capacities[idx] += 1; capacities[-1] -= 1

    return dict(capacities=capacities, penalties=penalties, F_rel_prob=F_rel_prob)   

Optimizing Bed Allocation

In [ ]:
# Run gradient descent for different bed totals (Takes a few minutes) 
# out165 = gradient_descent([55, 40, 30, 20, 20, 0])
# capacities165 = out165['capacities']
# out170 = gradient_descent([55, 40, 30, 20, 20, 5])
# capacities170 = out170['capacities']
# out180 = gradient_descent([55, 40, 30, 20, 20, 15])
# capacities180 = out180['capacities']

# Found capacities
capacities165 = [49, 28, 22, 16, 16, 34]
capacities170 = [48, 31, 25, 18, 14, 34]
capacities180 = [52, 35, 25, 16, 18, 34]
In [ ]:
out165 = {'capacities': [50, 32, 25, 14, 10, 34],
 'penalties': np.array([ 994.7 ,  992.51, 1023.78, 1016.06, 1031.28, 1054.56, 1037.07,
        1060.79, 1059.14, 1096.48, 1074.49, 1084.01, 1090.24, 1089.09,
        1106.09, 1110.83, 1113.56, 1110.67, 1121.34, 1120.74, 1126.26,
        1139.59, 1153.75, 1174.76, 1190.52, 1155.77, 1203.88, 1227.05,
        1269.71, 1258.28, 1284.74, 1305.39, 1327.39, 1327.1 , 1349.44]),
 'F_rel_prob': np.array([1.        , 0.96497114, 0.93127537, 0.89639315, 0.8651794 ,
        0.83104137, 0.7953772 , 0.76118428, 0.72714368, 0.69326257,
        0.66094461, 0.62798531, 0.59402437, 0.5632216 , 0.53226308,
        0.49888117, 0.46506489, 0.43620332, 0.40787952, 0.37578122,
        0.34587529, 0.31573004, 0.2874799 , 0.2606288 , 0.23393766,
        0.20971284, 0.18552666, 0.16272165, 0.14025691, 0.11959395,
        0.09943001, 0.0831262 , 0.06695472, 0.05308433, 0.04153909])}

Theoretical Optimal Number of Beds in F

In [ ]:
def erlangB_formula(m):
    if m == 0:
        return 1

    lam = 13; s = 2.2
    A = lam*s
    def power_div_factorial(i):
        return np.exp(i*np.log(A) - np.sum(np.log(range(1, i+1))))
    return power_div_factorial(m)/(np.sum([power_div_factorial(i) for i in range(m+1)]))
In [ ]:
fig = go.Figure()
fig.add_trace(go.Scatter(x=np.arange(35), y=[erlangB_formula(m) for m in range(35)], mode='lines', name='Erlang B'))
fig.add_trace(go.Scatter(x=np.arange(35), y=out165['F_rel_prob'], mode='lines', name='Simulation'))
fig.update_layout(title='Probability of relocating in ward F', xaxis_title='Number of beds in ward F', yaxis_title='Probability', width=800)
fig.show()

fig = go.Figure()
fig.add_trace(go.Scatter(x=np.arange(35), y=out165['penalties'], mode='lines', name='Simulation'))
fig.update_layout(title='Penalty as a function of number of beds in ward F', xaxis_title='Number of beds in ward F', yaxis_title='Penalty', width=800)
fig.show()

Estimate Rates for Admission, Relocation and Rejection

In [ ]:
def estimate_rates(samples, capacities):
    num_samples = len(samples)
    rates = np.zeros((num_samples, 3, 6))
    for i in range(num_samples):
        patients, events = samples[i]
        sim_out = simulate(patients, events, capacities)
        counts = sim_out['counts']
        rates[i, 0] = counts[ADMITTED]
        rates[i, 1] = counts[RELOCATED]
        rates[i, 2] = counts[REJECTED]

    rates /= rates.sum(axis=(1, 2))[:, None, None]
    rates = rates.mean(axis=0)
    rates = pd.DataFrame(rates.T, index=WARDS['names'], columns=['ADMITTED', 'RELOCATED', 'REJECTED'])
    return rates
In [ ]:
num_samples = 100
sim_setup_samples = [simulation_setup() for _ in range(num_samples)]

for n, cap in zip([165, 170, 180], [capacities165, capacities170, capacities180]):
    rates = estimate_rates(sim_setup_samples, cap)
    print(f'{n} beds:')
    print(rates.to_string(), '\n')
165 beds:
   ADMITTED  RELOCATED  REJECTED
A  0.206927   0.016383  0.029591
B  0.101809   0.043963  0.041178
C  0.056730   0.051762  0.029614
D  0.085913   0.014185  0.010681
E  0.039231   0.022740  0.022262
F  0.216036   0.006077  0.004917 

170 beds:
   ADMITTED  RELOCATED  REJECTED
A  0.208288   0.015399  0.029214
B  0.114944   0.040273  0.031734
C  0.065890   0.045475  0.026742
D  0.097142   0.008047  0.005590
E  0.036114   0.026882  0.021237
F  0.216036   0.006548  0.004446 

180 beds:
   ADMITTED  RELOCATED  REJECTED
A  0.225689   0.011821  0.015391
B  0.127633   0.034712  0.024605
C  0.070692   0.048506  0.018909
D  0.094211   0.010523  0.006044
E  0.049021   0.020295  0.014916
F  0.216036   0.007356  0.003638